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ABSTRACT 

We address the question of whether or not assembly bias arises in the absence of 
highly non-linear effects such as tidal stripping of halos near larger mass concentra- 
tions. Therefore, we use a simplified dynamical scheme where these effects are not 
modeled. We choose the punctuated Zel'dovich (PZ) approximation, which prevents 
orbit mixing by coalescing particles coming within a critical distance of each other. 
A numerical implementation of this approximation is fast, allowing us to run a large 
number of simulations to study assembly bias. We measure an assembly bias from 
60 PZ simulations, each with 512'^ cold particles in a 128ft-~^Mpc cubic box. The as- 
sembly bias estimated from the correlation functions at separations of ^ 5/i~^Mpc 
for objects (halos) at z = is comparable to that obtained in full N-body simula- 
tions. For masses 4 x lO^^/i^^M© the "oldest" 10% haloes are 3-5 times more cor- 
related than the "youngest" 10%. The bias weakens with increasing mass, also in 
agreement with full N-body simulations. We find that halo ages are correlated with 
the dimensionality of the surrounding linear structures as measured by the parameter 
(Ai -I- A2 -I- \?,)/{\\ + X2 + A|)^/^ where Xt are proportional to the eigenvalues of the 
velocity deformation tensor. Our results suggest that assembly bias may already be 
encoded in the early stages of the evolution. 
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1 INTRODUCTION 



A fundamental question in cosmology is the relation between 
the galaxy distribution and the underlying density field of 
the gravitationally dominant dark matter. According to the 
standard paradigm of structure formation, galaxies are har- 
bored in stable virialized objects (halos) made of dark mat- 
ter particles. An assumption that has been often made is 
that the clustering properties of halos depend on their mass 
alone. Although the assumption seems over-simplistic given 
the complexity of the hierarchical process of halo forma- 
tion, it is sustained by the excursion set theory (Bond et 
al. 1991; Lacey & Cole 1993; Mo & White 1996) and by re- 
sults of N-simulations of intermediate resolution (Lemson & 
Kauffmann 1999; Percival et al. 2003). Only recently Gao, 
Springel & White (2005) used a simulation of exceptionally 
large dynamical range (the Millennium Simulation, Springel 
et al. 2005) to show that the clustering of halos depend also 
on the their age, which is defined as the time since a halo 
acquired half of its current mass. They have found, that the 
"oldest" 10% of the halos with mass IQ^^ h^^ Mq are more 
than 5 times more correlated then the "youngest" 10% halos 
of the same mass. This assembly bias has been confirmed by 
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Harker et al. (2006) using marked correlation functions on 
the same simulation, and by Wechsler et al. (2006) and Jing 
& Mo (2006) using independent simulations. Wetzel et al. 
(2007) also found dependence of clustering on halo history, 
but only when using a different definition for the assembly 
redshift. 

We still lack a completely satisfactory explanation for 
the origin of assembly bias. For gaussian initial conditions, 
simple arguments based on the spherical collapse model ap- 
plied to narrow and broad initial density peaks which would 
collapse to halos of the same mass at the present time do 
predict an assembly bias, but with younger halos being more 
clustered than older ones. This trend of the bias is opposite 
to what is seen in simulations. Tidal stripping has been also 
invoked as a possible mechanism (e.g. Diemand, Kuhlen & 
Madau 2007). Because of mass stripping by the tidal gravi- 
tational field of the large mass concentration, nearby halos 
would have been of higher mass in a different environment. 
Therefore, these halos would have earlier formation times 
and would be more biased than halos of the same mass in 
the field. Avila- Reese et al. (2005) suggested tidal stripping 
as a mechanism responsible for generation of assembly bias 
in the high density regions whereas in low-density regions 
the cosmological initial conditions play a more important 
role. Maulbetsch et al. (2007) and Wang, Mo & Jing (2006) 



2 J. A. Keselman & A. Nusser 




16 32 48 16 32 48 64 

X [h-'Mpc] X [h-'Mpc] 



Figure 1. A visual comparison between results of the PZ approximation (right) and a PM N-body code (left) run with the same initial 
conditions for 256'^ particles in a box of 64Mpch~^ on the side. Top: particle (object for PZ) distributions in a slice of 3.2Mpch~^ in 
thickness (for the PM only a random subset of all the particles is shown). For the PZ, each object is represented as a filled circle with radius 
proportional to the mass. Bottom: contours maps of logj^g (density) in the same slice. The density fields are smoothed with a Gaussian 
window of a width of 1.125Mpc/i~^. The thin solid and dashed lines are density contours above and below the mean, respectively. Thick 
solid lines indicate mean density contours. The contour spacing is 0.18. 



and Desjacques (2007) suggest that the halo mass-accretion 
is less efficient in denser regions due to large scale tidal fields. 
However, the extent of this effect is difficult to assess. Here 
we examine whether the bias can, at least partially, arise 
in the quasi-linear evolution (i.e. over scales where the fiow 
is still laminar). In order to eliminate highly non-linear ef- 
fects such as tidal stripping, we adopt approximate meth- 
ods based on the Zel'dovich approximation (Zel'dovich 1970) 
where particles move on straight lines independent of the 
motion of other particles. The Zel'dovich approximation is 
an analytic solution to planar cosmological perturbations 
up to the stage where multi-streaming appears. For three 
dimensional perturbations, the approximation is a reason- 
able description of quasi-linear dynamics away from multi- 
streaming regions (e.g. Nusser et al. 1991). In order to ex- 
tend the applicability of this approximation beyond multi- 
streaming, we adopt the following scheme. Particles initially 
move in straight lines according to Zel'dovich, but they are 
merged together when they come within a critical distance 
of each other. This merging (sticking) produces an object 
with mass and linear momentum equal the total of its com- 
ponents. The critical distance is taken to depend on time 
like a diffusion length, as inspired by the adhesion approx- 
imation. This is known as the punctuated Zel'dovich (PZ) 
approximation (Fontana et al. 1995). This approximation 
is ideal for our purposes as it does not incorporate highly 
nonlinear effects and also it is fast and easy to implement. 



Further, it readily provides merging trees for individual ob- 
jects. 

The outline of the paper is as follows. In Sj2]we discuss 
the Zel'dovich approximation and our implementation of the 
punctuated Zel'dovich (PZ) scheme to describe the evolution 
in the quasi-linear regime. In |3]we present the PZ simula- 
tions and compare them with results from full dynamics. 
Our main results for halo biasing in the PZ simulations are 
presented in U We conclude with a general discussion in !j5] 



2 THE APPROXIMATE DYNAMICS 

According to the Zel'dovich approximation, the time depen- 
dent position, X, of a particle with initial (Lagrangian) co- 
ordinate q is 

X = q-f7?(t)0(q) , (1) 

where we work with comoving coordinates, the function D{t) 
is the gravitational growth rate of linear density fiuctuations, 
and 6{q) is a vector field which is a function of q only and 
is assumed to be derived from a potential. The physical pe- 
culiar velocity is v = a(f)dx/dt = aDO, where a{t) is the 
expansion factor of the cosmological background. The rela- 
tion (m yields a reasonable description of the gravitational 
evolution in the quasi-linear regime (e.g. Nusser et al. 1991; 
Weinberg & Gunn 1990), but fails near collapsed regions 
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where muhi-streaming occurs. The punctuated Zel'dovich 
(PZ) approximation is an extension of ([l}, which prevents 
the coasting away of particles in collapsed regions and yet 
preserves the simplicity of the Zel'dovich kinematics. In the 
PZ, one starts with equal mass particles located at a uni- 
form cubic grid at some initial time ti ^ 0. A particle is 
displaced from an initial position, q, at ti, according to ([1]) 
until it comes within a distance of dc of another particle. The 
two particles are then merged together to form an "object" 
having their total mass and linear momentum, and placed 
at their center of mass. The newly formed object is then 
moved according to the Zel'dovich relation ([ij with G deter- 
mined from momentum conservation. The scheme is applied 
to describe the further merging and evolution of objects (and 
particles) . 

We interpret the PZ in the framework of the adhe- 
sion approximation according to which (e.g. Shandarin 1991, 
Nusser & Dekel 1990) 



The viscous term modifies the Zel'dovich ansatz G = 
constant in regions with high velocity gradients and pre- 
vents orbit crossing. Viscosity affects the flow over (comov- 
ing) scales ^ V Du. Above these scales the flow is described 
by the usual Zel'dovich approximation. For dc oc i/D{t), 
the PZ is reminiscent of the adhesion approximation except 
that it ignores the details of the flow on scale ^ V Dv by 
coalescing objects which are within a distance dc of each 
other. Here we work with dc = V vD, as motivated by the 
adhesion approximation. 



3 THE SIMULATIONS 

We have run 60 PZ simulations, each with 512'' equal mass 
particles in a 128/i~^Mpc cubic box on the side. The initial 
conditions correspond to a random Gaussian realization of 
the Cold Dark Matter scenario in a flat universe without a 
cosmological constant. The dependence on the background 
density parameters comes through the initial power spec- 
trum but the dynamics is nearly independent of these pa- 
rameters when the linear growth factor is used as the time 
variable (Nusser & Colberg 1998). Therefore, apart from the 
effect of the initial power spectrum, the final result should be 
independent of the cosmological background. Thus the par- 
ticle mass is 4.3 x lQ^h~^ Mq. The dimensionless value of the 
Hubble constant is h = 0.73 and the rms value of the initial 
mass fluctuation in a sphere of 8h~^Mpc, as extrapolated to 
current time, is = 0.8. The particles are moved forward 
in time from z — 1000 to z = according to the PZ approx- 
imation with dc — V vD with u = lh~^Mpc^ which sets a 
spatial resolution of lh~^Mpc at the present time (D = 1). 
To further improve performance of the PZ, we smooth the 
initial velocity and density fields with a Gaussian window of 
width equal to a = l/i~^Mpc. 

Our time variable is the growth factor D and the time- 
step, dD, is such that |5max|d_D = O.ldc, where Smax is the 
speed of the fastest object. 

For purposes of history tracking, each object is assigned 
a unique ID. When objects merge, the newly formed objects 
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Figure 2. Density correlation functions computed from the PZ 
and PM simulations, as indicated in the figure. The mean of cor- 
relation functions from 60 PZ runs (each of 512^ particles in a 
box of 128 h~^Mpc on the side) is shown as the solid line. The 
attached error-bars are Icr. For this comparison we use a PM 
simulation of 256^ particles in a box size of 512h~^Mpc on the 
side. 

inherits the ID of the most massive progenitor. Object his- 
tories are tabulated in time slices separated by SD = 1/100. 

At the final time, the average number of objects per 
simulation is 7 x 10^ with 10* being more massive than 
4.3 X lO'^^h~^M0 which is the minimal halo mass we con- 
sider for the study of assembly bias. The number of halos in 
all 60 simulations is comparable to that in the millennium 
simulation. 

The PZ approximation is neither expected nor intended 
to model highly non-linear dynamics. Indeed, it is because it 
misses highly non-linear effects that we use it in this study. 
However, it is prudent to make a general comparison of our 
implementation of the PZ scheme with results from full dy- 
namics. To make a direct comparison, we have run the PZ 
scheme and a Particle-Mesh (Bertschinger & Gelb 1991) N- 
body code on the same initial conditions for 256^ particles 
in cubic box of 64h~^Mpc on the side. The initial conditions 
correspond to a flat CDM universe without a cosmological 
constant and erg — 0.8. Figure ([l| offers a visual impression 
of the difference between the final results from the PZ (pan- 
els to the right) and PM simulations (to the left) . In the top 
panels the final distribution of objects in the PZ approxi- 
mation is seen to follow closely the particle distribution in 
the PM code. This impression is further confirmed by the 
contour maps of the density fields in the bottom panels. The 
general agreement is impressive. The density fluctuations in 
the PZ are slightly of larger amplitude but this maybe due 
to cosmic variance. In figure ((2)| we also plot the density cor- 
relation functions obtained from PZ runs and a PM simula- 
tion. The solid line is the mean correlation computed from 
the mean of 60 PZ runs (each of 512'^ particles in a box of 
128h~^Mpc on the side), while the dashed is computed from 
the output of a single PM simulation of 256^ particles in a 
box of 512h-iMpc on the side. The la error-bars attached 
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to the PZ curve are estimated using the bootstrap method 
as follows. We generate 500 sets of simulations where each 
set contains 60 simulations picked randomly out of the 60 
original simulations (i.e. some of these simulations could be 
selected more than once) . For each set we compute the mean 
correlation and the errors are estimated as the standard de- 
viations between the mean correlations of the 500 sets. 

The bump in ^(r) in figure ([2]) at scales smaller than 
2h~^Mpc (also visible in figure [3]) is due to the finite resolu- 
tion in the simulations. Therefore, we will base our conclu- 
sion on correlations on scales larger than 2h~^Mpc 

Objects ("halos") in a PZ simulation are point-like and 
are identified using different criteria than halos in full N- 
body simulations. Therefore, we expect only a rough agree- 
ment between the mass functions of halos (number density 
versus mass) computed from PZ simulations and full dy- 
namics. We compared the abundance of objects versus mass 
in the PZ runs with the analytic predictions of Sheth & Tor- 
men (2002) and Press & Schechter (1974) for the halo mass 
function. The transfer function used in these predictions is 
taken from Bardeen et al. (1986) with a slope of n = 1 for 
the primordial power spectrum. 

Overall, there is only a qualitative agreement be- 
tween PZ and the analytic expressions. For masses 4.3 x 
lO"/i~^M0 the PZ simulation agrees with PS and ST. How- 
ever, for more massive haloes, the PZ overestimates abun- 
dance up to a factor of two for masses lQ^^h~^ Mq. The 
difference is reduced as we go to higher masses until it dis- 
appears at 6 X 10^'^/i~^Mq. At higher masses, PZ falls short 
of the analytic expressions by a factor which increases with 
mass. 



4 RESULTS 

The merging history of an object ("halo") in our imple- 
mentation of the PZ approximation is readily provided. 
We consider only halos containing more than 100 particles 
(4.3 X lO"/i"^M0) at the final time {z = 0) and define 
the formation time of a halo as the redshift at which it 
has acquired half of its final mass (Gao, Springel & White 
2005). We use the correlation functions to probe the clus- 
tering properties of halos. In figure|3]we plot the correlation 
functions, £,{x), as a function of separation, x, for halos in 
three mass ranges in the left, middle and right columns, re- 
spectively. The dashed (dotted) lines in the top, middle and 
bottom panels, respectively, correspond to 10%, 20% and 
30% oldest (youngest) halos. The solid lines in all panels are 
identical and represent the correlation function of the mass 
density field. In each panel the halo correlations are shown 
by two curves representing ±a deviations computed using 
the bootstrap method, as outlined in iJS] The dependence of 
the correlation function on the formation time is clear for 
all mass ranges shown in the figure. The bias persists even 
between the 30% youngest and 30% oldest halos. 

We use the difference between the correlation functions 
of old and young halos to quantify the assembly bias at 
various separations. We determine the bias parameter b in 
separation range {x,x + Aa;) by minimizing the quantity 
(Gao & White 2007) 



Figure (|4| shows the bias as a function of halo mass, for var- 
ious separations. For masses ^ 2 — 3 x lO^^M©, the bias is 
about 1.7 and is similar for all separations considered here. 
The error-bars are large at separations > 10h~^Mpc and we 
cannot detect an increase in the bias as claimed by Gao, 
Springel & White (2005). The bias weakens with increas- 
ing halo mass, but remains statistically significant only for 
the 10% old/young halos, at separations ^ 8h~^Mpc. The 
figure shows that the mass scale 2 — 3 x Mq marks 
a mass threshold above which assembly bias weakens, for 
all separations. This threshold is close to the non-linear 
mass scale Af* defined as the mass scale over which the 
rms of density fiuctuations is 1.69. For our initial conditions 
Af* « 5 X 10^^/i~U/q. 

Assembly bias may be caused by different environments 
of old and young halos. We have experimented with cross 
correlating the bias with several statistical measure of the 
environment. The most relevant measure that we find is the 
"dimensionality" of the density field in regions near halo 
particles at the initial time. This parameter is an indicator 
of the geometry of the structure developing at later times 
in those regions. We show here that halo ages are strongly 
correlated with the "dimensionality" of initial fluctuation 
field as defined by 



Al + A2 + A3 
VAfTAfTA 



(3) 



/ da; [log ^oid - log (6^ x Cyoung)]' 



where A; are the eigenvalues of the tensor dOi/dqm- At the 
centers of spherical, cylindrical and planar perturbations rj 
obtains the values 77 = \/3, \/2 and ?7 = 1, respectively. We 
have computed the mean value 7; as a function of distance 
from particles making up young and old halos. The results 
are plotted in figure Fig. (5) Solid and dotted lines, respec- 
tively, show rj for old and young halos (two lines representing 
zba deviations from the mean, calculated as explained in §3). 
This figure shows clearly that young haloes have an average 
higher dimensionality than old ones. 



5 CONCLUSIONS 

We have shown that assembly bias of halos persists even 
in a simplified description of gravitational dynamics like the 
punctuated Zel'dovich (PZ) approximation. The PZ approx- 
imation prevents the coasting away of particles in multi- 
streaming regions by coalescing objects that have come 
within a critical distance of each other. The PZ is fast, 
simple to implement, and readily provides object merging 
trees. This allows us to study assembly bias in a large num- 
ber of simulations (60 simulations, each of 512^ particles in 
a (128/i~^Mpc)'' cubic box). The magnitude of the bias is 
comparable to that found in full simulations. This implies 
that highly non-linear effects such as mass loss from halos 
in the vicinity of larger mass concentrations, may not be the 
dominant driver for assembly bias. 

We intend to apply the PZ scheme to the initial con- 
ditions used in the millennium simulation (Springel et al. 
2005) and compare the associated assembly bias with the 
resuh of Gao, Springel & White (2005). This will yield a bet- 
ter quantitative assessment of the role of highly non-linear 
effects. 

We have found a strong correlation between halo ages 
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Figure 3. Correlation functions for 10%, 20% and 30% oldest /youngest haloes. In each panel, the correlations are represented by curves 
corresponding to ztcr deviations. Dashed and dotted lines are for old and young halos, respectively. For comparison, the solid line in each 
panel shows the correlation function of the underlying mass- density. 



and the dimensionality of the nearby initial configuration. 
Young halos tend to form in regions of higher initial dimen- 
sionality than old halos. This is explained by the dependence 
of collapse time on dimensionality- a spherical perturbation 
collapses slower than a planar perturbation with the same 
initial density (Bertschinger & Jain 1994). Therefore, assem- 
bly bias would be explained if one could show that regions 
of lower dimensionality are more correlated than those of 
higher dimensionality. This problem could easily be tackled 
numerically using realizations of random gaussian fields. 
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Figure 5. The mean dimensionality parameter, ry, as a function of separation from particles making up young (dotted blue lines 
representing ±(t deviations from the mean) and old haloes (solid black lines), all for 10% oldest/youngest haloes. 



